
#include "Lattice.h"
#include <string>
#include <sstream>

using namespace std;

#define Nx			32
#define Ny			32

#define Rho			1.0
#define Nu			0.002
#define Se			1.0
#define Sepsilon	1.4
#define Sq			1.7
#define Snu			(2.0/(6.0*Nu+1))

#define U0			0.05
#define Kx			(2.0*M_PI/Nx)
#define Ky			(2.0*M_PI/Ny)
#define P0			0.0

#define Ntimes		60
#define Nglobals	1

int main(int argc, char *argv[]) {
	Lattice lattice;
	lattice.init(Nx, Ny, Rho, Se, Sepsilon, Sq, Snu);
	cout.precision(10);

	for(int i=0; i < Nx; i++)
		for(int j=0; j < Ny; j++) {
			long double u0x = -U0*cos(Kx*i)*sin(Ky*j);
			long double u0y = (Kx/Ky)*U0*cos(Ky*j)*sin(Kx*i);
// 			long double P = (-0.5)*U0*U0*cos(Kx*(i+j))*cos(Ky*(i-j));
// 			if (P < 0) cout << P << " " << i << " " << j << endl;
// 			long double _rho = 3.0*P;
			lattice.setU(i, j, u0x, u0y, Rho);
		}

// 	lattice.relaxRho();
	for(int n=0; n<Ntimes; n++) {
		if (n % int(Nglobals) == 0)
			cout << n << "\t" << lattice.energy()*4.0/(Nx*Ny*U0*U0) << endl;
		lattice.update();
	}
	return 0;
}
